arXiv: 1507.06348v3 [hep-th] 12 Jan 2016 


NRCPS-HE-24-2015 


Anosov C-systems 
and 

Random Number Generators 


George Savvidy 

Institute of Nuclear and Particle Physics 
Demokritos National Research Center, Ag. Paraskevi, Athens, Greece 


Abstract 

We are developing further our earlier suggestion to use hyperbolic Anosov C-systems 
for the Monte-Carlo simulations in high energy particle physics. The hyperbolic dynamical 
systems have homogeneous instability of all trajectories and as such they have mixing of all 
orders, countable Lebesgue spectrum and positive Kolmogorov entropy. These extraordi¬ 
nary ergodic properties follow from the C-condition introduced by Anosov. The C-condition 
defines a rich class of dynamical systems which span an open set in the space of all dynam¬ 
ical systems. The important property of C-systems is that they have a countable set of 
everywhere dense periodic trajectories and that their density exponentially increases with 
entropy. Of special interest are C-systems that are defined on a high dimensional torus. 
The C-systems on a toms are perfect candidates to be used for Monte-Carlo simulations. 
Recently an efficient algorithm was found, which allows very fast generation of long tra¬ 
jectories of the C-systems. These trajectories have high quality statistical properties and 
we are suggesting to use them for the QCD lattice simulations and at high energy particle 
physics. 

Dedicated to the memory of Professor Dmitri Anosov 



1 Introduction 


In the fundamental work on geodesic flows on closed Riemannian manifolds V n of negative 
curvature [I] Dmitri Anosov pointed out that the basic property of the geodesic flow on 
such manifolds is a uniform instability of all trajectories , which in physical terms means 
that in the neighbourhood of every fixed trajectory the trajectories behave similarly to the 
trajectories in the neighbourhood of a saddle point (see Fig. [I]). In other words, the 
hyperbolic instability of the dynamical system {T*} which is defined by the equations^] 

w = f(w ) (1-1) 


takes place for all solutions 6w = oo of the deviation equation 



ou 

w(t)=T t w 


( 1 - 2 ) 


in the neighbourhood of each phase trajectory w(t) = T f w, where w G IF™. 

The exponential instability of geodesics on Riemannian manifolds of constant negative 
curvature has been studied by many authors, beginning with Lobachevsky and Hadamard 
and especially by Hedlund and Hopf 121 E]. The concept of exponential instability of 
every trajectory of a dynamical system appears to be extremely rich and Anosov suggested 
to elevate it into a fundamental property of a new class of dynamical systems which he 
called C-system^} The brilliant idea to consider dynamical systems which have local and 
homogeneous hyperbolic instability of all trajectories is appealing to the intuition and has 
very deep physical content. The richness of the concept is expressed by the fact that the 
C-systems occupy a nonzero volume in the space of dynamical systems US 

Anosov provided an extended list of C-systems jT] . Important examples of the C- 
systems are: i) the geodesic flow on the Riemannian manifolds of variable negative curvature 
and ii) C-cascades - the iterations of the hyperbolic automorphisms of tori. 

Already from the intuitive definition of the C-systems it is clear that the C-systems 
have very strong instability of their trajectories and in fact the instability is as strong as 
it can be in principle [If. The distance between infinitesimally close trajectories increases 


1 It is understood that the phase space manifold W m is equipped by t he in variant Liouville measure [I]. 
2 The letter C is used because these systems fulfil the ”C condition” (1.3) [I]. 

3 This is in a contrast with the integrable systems where only a part of the trajectories remain on the 
invariant tori in accordance with the KAM theorem. 
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Figure 1: The integral curves in the case of saddle point x = —x, y — y. The t axis is 
the intersection of two planes of integral curves which are approaching it at t —>■ Too, the 
plane (x,t) and at t —¥ —oo, the plane (y,t). The rest of the integral curves are expanding 
at t —> Too, as well as at t —>■ — oo. A similar behaviour takes place in the neighbourhood 
of almost all trajectories of the C-system [T]. 


exponentially and on a compact phase space of the dynamical system this leads to the uni¬ 
form distribution of almost all trajectories over the whole phase space. Thus the dynamical 
systems which fulfil the C-condition have very extended and rich ergodic properties [T] . As 
such they have mixing of all orders, countable Lebesgue spectrum and positive entropy [I] . 

These properties follow from the theorems proven by Anosov and stating that the C- 
systems, with one exception]^ are in fact also Kolmogorov K-systems ME! . Thus the 
C-condition, in most of the cases, is a sufficient condition for the dynamical system to be 
a K-system as well. In this sense the C-systems provide extended and rich list of concrete 
examples of K-systems [1]. At the same time there are examples of K-systems which are not 
C-systems, therefore in the space of vector fields {f(w)} defining the dynamical systems 
w = f(w ), these two classes of dynamical systems are overlapping but are not identical 
(see also Footnote (§)■ 

The other important property of the C-systems is that in ’’between” the uniformly 
distributed trajectories there is a countable set of periodic trajectories. The set of points 
on the periodic trajectories is everywhere dense in the phase space manifold W m [T] . The 
periodic trajectories and uniformly distributed trajectories are filling out the phase space 
of a C-system in a way very similar to the rational and irrational numbers on the real line. 

Let us define the C-condition for the dynamical systems with discrete time jT|. A 
cascade on the m-dimensional compact phase space W m is induced by the diffeomorphisms 
T : W m — > IF" 1 . The iterations are defined by a repeated action of the operator {T n , —oo < 

4 These are the C-flows produced by some standard construction from the corresponding C-cascades [T] 
(see also Appendix). 
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a) 


b) 


Figure 2: The tangent vector oj G R w at the point w G IF 2 is decomposable into the sum 
R U' Xyj © Y w where the spaces X w and Y w are defined by the corresponding eigenvectors 
of the 2x2 matrix T (8.53). The automorphisms T induces the mapping of the tangent 
spaces T X w X'pyj . TY W Ytvj . 


It is contracting the distances on X w and expanding the 


distances on Y.„ 


n < +oo}, where n is an integer number. The tangent space at the point w G W m is 
denoted by R™ and the tangent vector bundle by lZ{W m ). The diffeomorphism {T n } 
induces the mapping of the tangent spaces T n : R™ —> R r f„ w . The C-condition requires 
that the tangent space R™ at each point w of the m-dimensional phase space W m of the 
dynamical system {T n } should be decomposable into a direct sum of the two linear spaces 
X k and Yf with the following properties [T]: 

Cl. K = (1.3) 

The dynamical system {T n } is such that : 

C2. a)\f n £\ < a\Z\e- cn for n > 0; \T n f\> b\f\e~ cn forn< 0, £ € X*, 
b)\T n r]\ > b\r]\e cn for n > 0; \T n r]\ < a\rj\e cn for n < 0, y G Y^, 

where the constants a,b and c are positive and are the same for all w G W m and all f G X*, 
y eY^. The length |...| of the tangent vectors £ and y is defined by the Riemannian metric 
ds on W m . 

The linear spaces X h w and Y l w are invariant with respect to the derivative mapping 
T n X l f = X!fn w , T n Y^ = Y^n w and represent the contracting and expanding linear spaces 
(see Fig§. The C-condition describes the behaviour of all trajectories T n aj on the tangent 
vector bundle oj G R r f . Anosov proved that the vector spaces Xf and Yf are continuous 
functions of the coordinate w and that they are the target vector spaces to the foliations Y k 
and E* which are the surfaces transversal to the trajectories T n w on W m (see Fig. [ 3 ]). The 
contracting and expanding foliations Y k , and Y l w are invariant with respect to the cascade 
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Figure 3: At each point w of the C-system the tangent space A™ is decomposable into 
a direct sum of two linear spaces Y l w and X’f. The expanding and contracting geodesic 
flows are 7 + and 7 “. The expanding and contracting invariant foliations Y} w and are 
transversal to the geodesic flows and their corresponding tangent spaces are Yf and Xf. 


T n in the sense that, under the action of these transformations a foliation transforms into 
a foliation and, in general, into a different one jT] . 

The interest in Anosov hyperbolic C-systems is associated with the earlier attempts 
to understand relaxation phenomena of the hard scattering spheres and the foundation 
of the statistical mechanics by Krylov [Sj. The hyperbolic systems demonstrate strong 
statistical properties and can help with the understanding of the appearance of turbulence 
in fluid dynamics [9], non-integrability of the Yang-Mills classical mechanics mm , as 
well as dynamical properties of gravitating N-body systems and the relaxation phenomena 
of stars in galaxies mm- The present article is devoted to the application of C-systems 
defined on a torus for the Monte-Carlo simulations of physical systems [T3j . The random 
number generators based on C-systems demonstrated high quality statistical properties 
mm and are implemented into the ROOT and GEANT software used at CERN for the 
LHC experiments [HI [45l ©j. 

In the next sections we shall consider the properties of the C-systems on a torus. The 
interest in the behaviour of the C-systems on a torus is associated the with the fact that 
they can be used to generate high quality pseudorandom numbers for the Monte-Carlo sim¬ 
ulations of physical systems. In the second section we shall define a linear diffeomorphism 
on a torus which fulfils the C-condition and shall describe the distribution properties of 
their periodic trajectories. In the third section we shall define the Kolmogorov entropy of a 
dynamical system and calculate its value for the diffeomorphisms on a torus. In the fourth 
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section we shall demonstrate that the integrals over the phase space can be calculated as 
a sum over points of the infinitely long periodic trajectories. In the fifth section we shall 
present recent progress in the efficient implementation of a particular C-system to generate 
pseudorandom numbers for Monte-Carlo simulations and their usage in high energy physics 

[32113311331135113B] . 

2 Torus Automorphisms 


Let us consider a dynamical system T with the phase space IT™ which is the m-dimensional 
torus appearing at factorisation of the Euclidean space E m with coordinates (uq, .... w m ) 
over an integer lattice. Then T can be thought of as a linear transformation of the Eu¬ 
clidean space E m which preserves the lattice L of points with integer coordinates. The 
automorphisms of the torus are generated by the linear transformation 

Wj —> ^ Tj jWj, (mod 1), (2.4) 


where the integer matrix T has a determinant equal to one DetT =1. In order for the 


automorphisms of the torus (2.4) to fulfil the C-condition it is necessary and sufficient that 
the matrix T has no eigenvalues on the unit circle. Thus the spectrum {A = Ai,..., A m } of 
the matrix T is chosen to fulfil the following two conditions [3]: 


1) DetT — A 1 A 2 _A m — 1 

2) |Ai| 7 ^ 1 , Vi (2.5) 


Because the determinant of the matrix T is equal to one, the Liouville’s measure d/v, = 
dx\...dx m is invariant under the action of T. The inverse matrix T~ l is also an integer 


matrix because DetT = 1. Therefore T is an automorphism of the torus W m onto itself. 


This automorphism has a fixed point p corresponding to the origin of E r 


The above conditions (2.5) on the eigenvalues of the matrix T are sufficient to prove 
that the system belongs to the class of Anosov C-system^] . For that let us divide the 


5 If within the eigenvalues of the matrix T there are no roots of unity then the automorphism T of a 
torus defines a K-system m ■ This means that some of the eigenvalues of T can be on a unit circle on 
the contrary to the C-conditions ( |2.5| ) . It follows then that the K-systems occupy a different region in the 
space of dynamical systems {w = f(w)}. In two and three dimensions the condition ( |2.5| ) is equivalent to 
the absence of eigenvalues equal to the root of unity, thus in these cases the K- and C- conditions coincide, 
but in dimensions larger that three (m > 4) there are K-systems which are not C-systems pQ. 
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eigenvalues of the matrix T into the two sets {A a } and {A^} with modulus smaller and 
larger than one: 


0 < |A„| < 1 < |A„|. 


( 2 , 6 ) 


Consider two family of planes X = {X a } and Y = {Yp} which are parallel to the cor¬ 
responding eigenvectors {e a } and {ep} . The derivative map T maps these planes into 
themselves contracting the points on the X a plane \ a times and expanding the points on 
the Yp plane A p times (see Fig. |4]): 


m < |A a | Id, 

\T1]\ > |A S | |>)|, r, e Y e . 


(2.7) 


It follows then that the automorphism T fulfils the C-condition (1.3) with a=b=l, e c = A, 
where 

A = min f - - min |A^| ) . (2.8) 

\max |A a | / 

In other words, the invariant planes of the matrix T , for which the eigenvalues are inside 
or outside of the unit circle define the contracting and expanding invariant spaces of the 


C2-condition (1.3), so that the phase trajectories of the dynamical system are expanding 


and contracting under the transformation T at the exponential rate. The foliations S a 
and Sp which are orthogonal to the phase trajectories are divided into contracting and 
expanding foliations and the tangent planes to the foliations comprise the contracting and 
expanding planes {X a } and {Yp}. The foliations can be constructed by filling out the 
Euclidean space E m by planes parallel to {X a } and {Yp} and then projecting them back 
to the torus W m . 

Let us now consider the image of the set F of W m under the action of the transformation 
T n . It will be a hyperbolic rotation of F: expanding its points in the ep directions in 
\Xp\ n = e nln l A/3 l rate and contracting in the e a directions in |A a | n = e -nln ^«T rate. Therefore 
at large n the image of the set T n F will be represented by a long stripe stretched in the 
direction of the largest eigenvalue A p and uniformly distributed over the phase space IT'" 
in a way similar to the uniformly distributed trajectory of the equation w = fl on a torus 
with non-resonant vector Lrl Thus, for large n, the images T n F are crossing an arbitrary 


6 The vector Q is non-resonant if the equation nXh + ... + n rn il m = 0 has no solutions with integer 
(ui, * - ., u m ). 
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Figure 4: The matrix T (8.53) has two irrational eigenvalues A a = (3 —\/5)/2 and Xp = (3+ 
V / 5)/2. The corresponding eigenvectors {e a } and { ep} define two families of parallel lines 
{X a } and {Yp} on the two-dimensional torus which are invariant under the automorphisms 
T. The automorphism T is contracting the distances between points on the lines belonging 
to the set {X a } and expanding the distances between points on the lines belonging to the 
set {Yp}. The a) depicts the parallel lines of the sets {X a } and {Yp} and b) depicts their 
positions after the action of the automorphism T. 


region G C W m proportional to its volume n(G): 


lim n(T n F n G) = n{F)n(G). (2.9) 

n—> oo 


This equation represents a mixing property of the dynamical system T. If one introduces a 


set of functions {g(w)} on the phase space W m , then one can represent the equation (2.9) 
in the form EZ0 


lim 

n—> oo 


f(T n w) g(w)d/j,(w) 


f{w)dn{w ) / g(w)dfj,(w). 


( 2 . 10 ) 


The above property is known in physics as the factorisation or decay of the two-point 
correlation functions. Because the C-systems have mixing of all orders [lj it follows that 
the C-systems are exhibiting the decay of the correlation functions of any order. 

A strong instability of trajectories of a dynamical C-system leads to the apperance of 
statistical properties in its behaviour int As a result of exponential instability of the 
motion on the phase space the time average of the function g(w) on W m 


1 

N 


N -1 

£ g(T"w) 

n =0 


( 2 . 11 ) 


behaves as a superposition of quantities which are statistically weakly dependent. Therefore 
for the C-systems on a torus it was demonstrated that the fluctuations of the time averages 


7 The equation (2.9) follows from (2.10) if one considers phase space functions which are equal to one 
on all points belonging to the sets F and G and to zero for points outside the sets. 
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(2.11) from the phase space integral 


iw " 


g(w)dw 


( 2 . 12 ) 


multiplied by s/N have at large IV —» oo the Gaussian distribution HZ]: 

i r 


lim a 

TV—>-00 


\w : W ( i ^ S (T”u,) - f 

{ V n=0 


g(v)dv < 2 > = 


\/2na 

q J — c 


e (2.13) 


The importance of the multiplication by the factor \[N can be understood as follows. 
The difference in the bracket has an upper bound in terms of the Kolmogorov discrepancy 
Dn(T) [151132j : 

D n (T) 


i iv_1 r 

- £ s(T"w) - / 9 (v)di 
7V n =0 Jwm 


< C 


N 


(2.14) 


where C is a constant and D^{T) grows as \f~N. Therefore after multiplication of the 
quantity in the bracket by \fN it is bound by a constant. The theorem states that the 


measure of points of W m which fulfil the inequality (2.13) have Gaussian distribution. 

Our next aim is to demonstrate that C-systems have a rich variety of periodic trajecto¬ 
ries, which essentially depends on the entropy of the C-systems Pl20l[2ll[23l[22[25l[26]. 
We shall provide below a short introduction to the definition of the Kolmogorov entropy 
and shall calculate its value for the different authomorphisms of a torus. 


3 The Entropy of the Automorphisms T 

A class of dynamical systems was introduced by Kolmogorov in 0E! which he called 
quasi-regular or K -systems and defined the notion of entropy for such systems. Consider 
dynamical systems with discrete time. Let a = {Ai} ieI ( I is finite or countable) be a 
measurable partition of the phase space W m (W m is equipped with a positive measure /.(), 
that is 

\ |J AP = 0, niAi P| Aj) = 0, i ± j , (3.15) 

iei 

then one can define the entropy of the partition a as 

h(a) = ~y^/r(As)ln//(Aj). (3.16) 

iei 
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It follows that if two partitions cq and « 2 differ by a set of zero measure then their entropies 
are equal. The refinement partition a 


a = ai V «2 V ... V ft*; 


(3.17) 


of the collection of partitions ai ,..., is defined as the intersection of all their composing 
sets Ap. 

« = U1 Ai | Ai G for all i}. (3.18) 

iei 

The entropy of the partition a with respect to the automorphisms T is defined as a limit 

nummmm-- 


, . h(a V Ta V ... V T n ~ l a) 

h(a, T) = lim —--- n = 1,2, 


n 


(3.19) 


This number is equal to the entropy of the refinement /?: 


f3 = a V Ta V ... V T n 1 o;, 


(3.20) 


which was generated during the iteration of the partition a by the automorphism T. Finally 
the entropy of the automorphism T is defined as a supremum: 


h(T) = sup h(a, T), (3.21) 

M 

where the supremum is taken over all finite measurable partitions {a} of W m . From the 
definition it follows that the calculation of this number for a given dynamical system seems 
extremely difficult. The theorem proven by Kolmogorov tells that if one finds the so 
called ’’generating partition” (3 for the automorphisms T which has the property 

OO 

U T n M((3) = M, (3.22) 

Tl— — OO 

where A4(/3) is the sigma-algebra of the partition f3 = {Bi} iG i and Ai is the sigma-algebra 
of all measurable sets of W m , then 


h(T) = h(/3,T). 


(3.23) 


8 A sigma-algebra M.(a) on W 171 is a family of subsets {A^ e / E a of the partition cp which is closed 
under union operations of countably many sets and complement 
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Which means that the supremum in (3.21) is actually reached on a generating partition (3. 
In some cases the construction of the generating partition (3 allows an explicit calculation 
of the entropy of a given dynamical system [22, 28j. 

There are also alternative ways of calculating the entropy of a C-system. The most 
convenient for us is the integration over the whole phase space of the logarithm of the 
volume expansion rate A (w) of a /-dimensional infinitesimal cube which is embedded into 
the foliation E[ (1 . The derivative map T maps the linear space Yj, into the Y^ w and if the 
rate of expansion of the volume of the /-dimensional cube is A(w), then [U 7, UHl 22, [2E] 

h(T) = ( In A (w)dw. (3.24) 

Jw m 


Here the volume of the W m is normalised to 1. For the automorphisms on a torus (2.4) 
the coefficient A (w) does not depends of the phase space coordinates w and is equal to the 


product of eigenvalues {A^} with modulus larger than one (2.6): 


am= n A*)- 

0=i 


(3.25) 


Thus for the Anosov automorphisms on a torus (2.4), (2.5) from (3.24) and (3.25) we can 
calculate the entropy, which became equal to the sum: 


h(T)= HA/sl- (3-26) 

|A/4>1 

It is interesting to know what will happen if the authomorphism of a torus is defined 
by the k’th power T k of the matrix T. In that case it will generate a cascade of the form 
{(T fc )"} and the question is what is the entropy of this authomorphism. The answer can 
be found by calculating the eigenvalues of the new matrix, therefore 


h(T k ) = ^ IA/ = k In \Xp\ = k h(T). (3.27) 

I-M>1 |Ag|>l 

As one can see, the entropy increases linearly with the power of the matrix k and, as 
we shall see in the next sections, this is a very ’’expensive” and time consuming way of 
increasing the entropy of the C-systems used for Monte-Carlo calculations. 

In the next section we shall present the alternative derivation of the above results, 
which will demonstrate the connection of the entropy with the variety and richness of the 
periodic trajectories of the C-systems [H [21 [251 HU- 
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4 Periodic Trajectories of Anosov C-systems 


Let us demonstrate now that the C-systems have a countable set of periodic trajectories 
[I]. The E m cover of the torus W m allows to translate every set of points on torus into a 
set of points on Euclidean space E m and the space of functions on torus into the periodic 
functions on E m . To every closed curve 7 on a torus corresponds a curve 0 : [0,1] —> E m 
for which 0 ( 0 ) = 0 ( 1 ) mod 1 and if 0 ( 1 ) — 0 ( 0 ) = (pi, then the corresponding 

winding numbers on a torus are p t G Z. 

All trajectories with rational coordinates (wi ,..., w m ), and only they, are periodic tra¬ 
jectories of the automorphisms of the torus ( |2.4 ). Let us fix the integer number N, then 
the points on a torus with the coordinates having a denominator N form a finite set 


{pi/N, ...,p m /N}. The automorphism (2.5) with integer entries transform this set of points 
into itself, therefore all these points belong to periodic trajectories. Let w = (wi, 
be a point of a trajectory with the period n > 1. Then 

T n w = w+p, (4.28) 

where p is an integer vector. The above equation with respect to w has nonzero determi¬ 
nant, therefore the components of w are rational. 

Thus the periodic trajectories of the period n of the automorphism T are given by the 


solution of the equation (4.28), where p G Z m is an integer vector and w = (uq, G 

W m . As p varies in Z m the solutions of the equation (4.28) determine a fundamental 
domain D n in the covering Euclidian space E m of the volume p(D n ) = l/\Det(T n — 1)|. 
Therefore the number of all points N n on the periodic trajectories of the period n is given 
by the corresponding inverse volume [ 20 , EH EH EH 25]: 

m 

(4.29) 


N n = \DetCr -1)1 = 1 JpA” - 1)|. 


2=1 


Using the theorem of Bowen |25l ES! which states that the entropy of the automorphism T 


can be represented in terms of N n defined in (4.29): 

h(T) = lim - In A^ n 

n—t 00 Jl 


(4.30) 


one can derive the formula for the entropy (3.26) for the automorphism T in terms of its 
eigenvalues: 

-I m 

h(T)= lim -ln(| JJ(A?- 1)|) = ^ ln|A^|. (4.31) 

n—¥00 n. A A L —' 

I a /3|>! 


2=1 
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Let us now define the number of periodic trajectories of the period n by i r(n). Then the 
number of all points N n on the periodic trajectories of the period n can be written in the 
following form: 

N n = Y MO, (4-32) 


l divi n 


where l divides n. Using again the Bowen result (4.30) one can get 

N n = Y MO ~ e nh(T) - (4.33) 

l divi n 

This result can be rephrased as a statement that the number of points on the periodic 
trajectories of the period n exponentially grows with the entropy. 

Excluding the periodic trajectories which divide n (for example T n w = T l ' 2 (T l> w). i 
where n = l-Jo and T li w = w ) one can get the number of periodic trajectories of period n 
which are not divisible. For that one should represent the nin) in the following form: 


M) = ^( Y MO- X! MO) 


(4.34) 


l divi n 


l divi n, l<n 


and from (4.34) and (4.33) it follows that 

a nh(T) y' 

divi n , l<n 


7 r (n) 


(1- 


l 7r(Z). 


a nh(T) 


(4.35) 


n ' Ei divi „ 1 n 

because the ratio in the bracket is strictly smaller than one. This result tells that a system 
with larger entropy Ah = h(Tj) — h(T 2 ) > 0 is more densely populated by the periodic 
trajectories of the same period n: 

7Ti(n) 


7r 2 (n) 


~ e 


n Ah 


(4.36) 


states that 


The next important result of the Bowen theorem 

f(w)dfj,(w) = lim -L Y 

v. -4oo /V._ * ^ 


I 


(4.37) 




where T n is a set of all points on the trajectories of period n. The total number of points 
in the set T n we defined earlier as N n . 

This result has important consequences for the calculation of the integrals on the man¬ 


ifold W m , because, as it follows from (4.37), the integration reduces to the summation 
over all points of periodic trajectories. It is appealing to consider periodic trajectories of 
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the period n which is a prime number. Because every infinite subsequence of convergent 


sequence converges to the same limit we can consider in (4.37) only terms with the prime 
periods. In that case N n = mr(n) and the above formula becomes: 


1 


7r(n) n—1 


/ f(w)dn(w ) = lim x / 

/ld/ m n^oo UTrln) *— 

J 7=1 i —I 




(4.38) 


where the summation is over all points of the trajectory T l Wj and over all distinct tra¬ 
jectories of period n which are enumerated by index j. The Wj is the initial point of the 
trajectory j [j From the above consideration it follows that the convergence is guaranteed 
if one sums over all trajectories of the same period n. One can conjecture that all 7r(n) 


trajectories at the very large period n contribute equally into the sum (4.38), therefore the 


integral (4.38) can be reduced to a sum over fixed trajectory 

n—1 


i£/(TV). 

2=0 


(4.39) 


This reduction seems plausible, but even in that case the calculation of (4.39) is impractical, 
because the periods of the C-systems which are used for Monte-Carlo calculations can be 
as large as n ~ lO 4000 jqjj The Monte-Carlo calculation of the integrals is therefore 
performed on some parts of the length N of the long trajectories [29, EDI ED E2j- In that 
case the approximation of the integrals is regulated by the Kolmogorov discrepancy Dj^{T) 


in the formula (2.14). The behaviour of D^{T) as a function of N crucially depends on 
the properties of the dynamical system T because the discrepancy grows slower for the C- 
systems periodic trajectories which are distributed uniformly and everywhere dense in the 
phase space W m [DE]. Therefore it was suggested that Anosov C-systems [lj, defined on a 
high dimensional torus, are excellent candidates for the pseudo-random number generators 
na. In the next section we shall demonstrate how these results can be used for practical 
calculations and especially in high energy particle physics. 


5 Anosov C-systems for Monte-Carlo Computations 

Modern powerful computers open a new era for the application of the Monte-Carlo Method 
[29, IM, 31, EH IE i36, 37;, 38j for the simulation of physical systems with many degrees 

9 It appears to be a difficult mathematical problem to decide whether two vectors w\ and W 2 belong to 
the same or to distinct trajectories. 
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of freedom and of higher complexity. The Monte-Carlo simulation is an important com¬ 
putational technique in many areas of natural sciences, and it has significant application 
in particle and nuclear physics, quantum physics, statistical physics, quantum chemistry, 
material science, among many other multidisciplinary applications. At the heart of the 
Monte-Carlo (MC) simulations are pseudo Random Number Generators (RNG). 

Usually pseudo random numbers are generated by deterministic recursive rules CSU291 
2EE ED E2J- Such rules produce pseudorandom numbers, and it is a great challenge to 
design pseudo random number generators that produce high quality sequences. Although 
numerous RNGs introduced in the last decades fulfil most of the requirements and are 
frequently used in simulations, each of them has some weak properties which influence the 
results [HS] and are less suitable for demanding MC simulations which are performed for the 
high energy experiments at CERN and other research centres. The RNGs are essentially 
used in ROOT and GEANT software for the design of the experiments and statistical 
analysis of the experimental data [33] , 

In order to fulfil these demanding requirements it is necessary to have a solid theoretical 
and mathematical background on which the RNG’s are based. RNG should have a long 
period, be statistically robust, efficient, portable and have a possibility to change and 
adjust the internal characteristics in order to make RNG suitable for concrete problems of 
high complexity. 

In [fBJ it was suggested that Anosov C-systems pQ, defined on a high dimensional torus, 
are excellent candidates for the pseudo-random number generators. The C-system chosen 
in was the one which realises linear automorphism T defined in (2.4) [^J 

N 

Wj —> TjjWj, (mod 1). (5.40) 

j= 1 


10 For convenience in this section the dimension m of the phase space W m is denoted by N. 
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A particular matrix chosen in [[23] was defined for all N > 3: 

(\ i i i i A 
1 2 
1 3 4" s 


T = 


1 4 


1 

1 

2 

3 


1 ...11 

1 ...11 

1 ...11 

2 ...11 


V 1 


N N- 1 N - 2 ... 3 


(5.41) 


2 / 


The matrix T of integer numbers is of the size N x N and has determinant equal to one. 
It is defined recursively, since the matrix of size N + 1 contains in it the matrix of the 
size N. The only variable entry in the matrix is T 32 = 3 + s, where s should be chosen 
such that to avoid eigenvalues lying on a unit circle. In order to generate pseudo-random 
vectors w n = T n w, one should choose the initial vector w = (uq...., w m ), called the “seed”, 
with at least one non-zero component to avoid fixed point of T, which is at the origin. 

The eigenvalues of the T matrix (5.41) are widely dispersed for all N, see Fig. [5] from 
reference m- The spectrum is ’’multi-scale”, with trajectories exhibiting exponential 
instabilities at different scales [T5]. The entropy of T satisfies the bound 


- N In 4 < h < N In 4, 
9 


(5.42) 


and tells that the entropy of T is increasing linearly with its size N. In the light of the 



Figure 5: The logariphm of the absolute value of the eigenvalues of the MIXMAX matrix 


for the size N=256. The area under the curve is proportional to the entropy (3.26). 


discussion in the previous section this means that there is a rich set of periodic trajectories 
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with large periods and their variety increases exponentially with the size N of the matrix 
T. For the matrix of the size N = 256 the actual value is h ~ 164.4 thus 

^164 n 


7 r (n) 


n 


(5.43) 


The Monte-Carlo simulations of a Markov chain should be done with a random number 
generator whose auto-correlation time To = 1/fi is much smaller than the auto-correlation 
time of the Markov chain 


In a typical computer implementation of the authomorphism (8.48) the initial vector 
will have rational components iq = a,/p, where a, and p are natural numbers. Therefore 
it is convenient to represent u, by its numerator a* in computer memory and define the 
iteration in terms of a* [Mj : 

N 

ai —y Tij aj mod p. (5.44) 

j= 1 

If the denominator p is taken to be a prime number [31], then the recursion is realised on 
extended Galois field GF\p N ] [39] 40J and allows to find the period of the trajectory n in 
terms of p and the properties of the characteristic polynomial P(x ) of the matrix T [34j. 
If the characteristic polynomial P(x) of some matrix T is primitive in the extended Galois 
field GF\p N ], then [34, 40, 41] 

p N — 1 


T q = Pq I where q = 


p — 1 


(5.45) 


where p 0 is a free term of the polynomial P(x) and is a primitive element of GF\p\. Since 
our matrix T has po = DetT = 1, the polynomial P(x) of T cannot be primitive. The 
solution suggested in [42] is to define the necessary and sufficient conditions for the period 
q to attain its maximum, they are: 


N _I 

1. T q = I (mod p ), where q = pprp 

2 . T q ' r 7 ^ I (mod p), for any r which is a prime divisor of q . 

The first condition is equivalent to the requirement that the characteristic polynomial is 
irreducible. The second condition can be checked if the integer factorisation of q is available 


, then the period of the sequence is equal to (5.45) and is independent of the seed. There 
are precisely p— 1 distinct trajectories which together fill up all states of the GF[p N ] lattice: 

q (p — 1) = p N — 1. (5.46) 
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In [32] the actual value of p was taken as p = 2 61 — 1, the largest Mersenne number that 
fits into an unsigned integer on current 64-bit computer architectures. For the matrix of 
the size N = 256 the period in that case is q ps lO 4600 . 

The algorithm which allows the efficient implementation of the generator in actual com¬ 
puter hardware, reducing the matrix multiplication to the O(N) operations was found in 
[42]. The other advantage of this implementation is that it allows to make ’’jumps” into 
any point on a periodic trajectory without calculating all previous coordinates on a trajec¬ 
tory, which typically has a very large period q ~ lO 4600 . This MIXMAX random number 
generator is currently made available in a portable implementation in the C language at 
hepforge.org [33] and was implemented into the ROOT and GEANT projects at CERN 
[3311351135] , 

6 Conclusion 

The author had the privilege to visit Professor Dimitry Anosov at the Steklov Institute in 
1986 to discuss and learn his fundamental work on hyperbolic dynamical systems. Pro¬ 
fessor Anosov’s impressive simplicity and clear explanation of complicated mathematical 
constructions was an unforgettable experience. 
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8 Note Added 


In [T] Anosov demonstrated how any C-cascade on a torus can be embedded into a certain 


C-flow. The embedding was defined by the identification (8.47) and the corresponding 


C-flow on a smooth Riemannian manifold W m+1 with the metric (8.55) was defined by 
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w -► Tw 



Figure 6: The identification of the W 2 x {0} with W 2 x {1} by the formula (w, 1) = (Tw, 0) 
of a cylinder W 2 x [0,1], where [0,1] = {u | 0 < u < 1}. The resulting compact manifold 
W 3 has a bundle structure with the base S 1 and fibres of the type W 2 . The manifold W 3 
has the local coordinates w = (w l ,w 2 ,u) . 


the equations (8.49). We are interested here to analyse the geodesic flow on the same 
Riemannian manifold W m+1 . The geodesic flow has different dynamics (8.56 ) and as we 


shall demonstrate below has very interesting hyperbolic components different from (8.49). 


Let us consider a C-cascade on a torus W m (defined in section two) and increase its 
dimension m by one unit constructing a cylinder W m x [0,1], where [0,1] = { u | 0 < u < 1}, 
and identifying W m x {0} with W m x {1} by the formula: 


(w, 1) = (Tw, 0). 


Here T is diffeomorphism (2.4): 


w —» 


^Tijwfl (mod 1). 


(8.47) 


(8.48) 


The resulting compact Riemannian manifold W m+1 has a bundle structure with the base 
S 1 and fibres of the type W rn . The manifold W m+1 has the local coordinates w = 
(w l ,..., w m , u ) . 

The C-flow T f on the manifold W m+1 is defined by the equations |[TJ 


dw 1 n dw m n du 1 
~dt = ~df = ’ dt = L 

1 


(8.49) 


For this flow the tangent space R™ +L can be represented as a direct sum of three subspaces 


- contracting and expanding linear spaces X^,Y^ and Z , ; ,: 

RT+ 1 = Xi © Yi © Z, r ,. 


(8.50) 
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The linear space X*~ is tangent to the fibre W m x u and is parallel to the eigenvectors 
corresponding to the eigenvalues which are lying inside the unit circle 0 < |A a | < 1 and 
is tangent to the fibre W m x u and is parallel to the eigenvectors corresponding to the 
eigenvalues which are lying outside of the unit circle 1 < |A ( g|. Zy, is collinear to the phase 


space velocity ( |8.49[ ). Under the derivative mapping of the ( |8.49[ ) the vectors ( |8.57[ ) from 
Xl and Y‘ are contracting and stretching: 


\f ‘ml = \\ |x'i|, ir'x^l = A| M- 


(8.51) 


This identification of contracting and expanding spaces proves that (8.49) indeed defines 
a C-flow [T]. 

It is also interesting to analyse the geodesic flow on a Riemannian manifold V n = W m+1 
(n=m+l). The equations for the geodesic flow on V n 


d 2 w p dw u dw p 

+ n —^- = o 


dt 2 


up dt dt 


(8.52) 


are different from the flow equations defined by the equations (8.49) and our goal is to 
learn if the geodesic flow has also the properties of the C-flow. The answer to this question 
is not obvious and requires investigation of the curvature structure of the manifold V n . If 
all sectional curvatures on V n are negative then geodesic flow defines a C-flow [T]. 

For simplicity we shall consider a two-dimensional case which is defined by the matrix 


T = 



(8.53) 


The metric on the corresponding manifold U 3 can be defined as ra 

ds 2 = e 2u [Xi dw 1 + (1 — Ai )dw 2 ] 2 + e 2u [A 2 dw 1 + (1 — A 2 )dw 2 ] 2 + du 2 = g^dw p dw l/ , 


where 0 < A 2 < 1 < Ai are eigenvalues of the matrix (18.53) and fulfil the relations 


AiA 2 = 1, Ai + A 2 = 3. The metric is invariant under the transformation 

w 1 = 2w' 1 — w' 2 , w 2 = —w'i + w' 2 , u — u—1 (8.54) 

and is therefore consistent with the identification ( |8.47| ). The metric tensor has the form 

X 2+2u + A2 +2 “ (1 - Ai)A} +2 “ + (1 - A 2 )A2 +2u o\ 




9^{u) = 


X 1 

i 1 —|— 2 u 


(1 - A!)AC 2u + (1 - A 2 )A^ +2u (1 - X 1 fX 2u + (1 - A 2 ) 2 A 2 “ 0 


V 


0 


0 


(8.55) 


V 
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and the corresponding geodesic equations take the following form: 


(Ai — lllnAi • (Ai — l)lnAi • 

w 1 + 2-^-——- -w l u - 4-^—-—- -w 2 u = 0 

Ai + 1 Ai + 1 

-2 0 (^i — 1) In Ai • (Ai —l)lnAi • 

W — 2 -:--- W 2 U — 4 - T -:- W L U = 0 


Ai + 1 


Aid 1 


(8.56) 


» (l-Ar 4 )lnA l m , ml ^ qi + A^KAi-lllnA^ ., 


il + 


y2u+2 


(l-Af +2 )(A 1 -l) 2 lnA 1 

y2u-\-2 


y2u+2 




W 2 W 2 = 0. 


One can become convinced that these equations are invariant under the transformation 


(8.54). In order to study a stability of the geodesic flow one has to compute the sectional 


curvatures. We shall choose the orthogonal frame in the directions of the linear spaces 
X1 ,, I)] and Zy,. The corresponding vectors are: 


Vi — (Ai — 1, Ai, 0), V 2 — (A 2 — 1, A 2 , 0), v$ — (0, 0,1) 


(8.57) 


and in the metric (8.55) they have the lengths: 


hi 2 = (Ai - A 2 ) 2 A 2 “, \v 2 \ 2 = (Ai - A 2 ) 2 Af, H 2 = 1. 


2 \ 2 u 


(8.58) 


The corresponding sectional curvatures are: 


v _ Rnv\pViV%V$V% 

R1 2 — - 1 -t -fy- 

\Vi A v 2 \ 2 

rv _ R ^\P V 1 V 3 V 1 V 3 _ 

R-13 — —|-t -ry— — 

\Vi A v 3 \ 2 

T S Ftpv\p V 2 V 3 V 2 V 3 

R 23 ~ -i-7-19- — 

\v 2 A u 3 | 2 


= In 2 Ai 
- In 2 A 2 


— In 2 Ai. 


(8.59) 


It follows for the above equations that the geodesic flow is exponentially unstable on the 
planes (1,3) and (2,3) and is stable in the plane (1,2). This behaviour is dual to the flow 
( 8.49[ ) which is unstable in (1,2) plane and is stable in (1,3) and (2,3) planes. The scalar 
curvature is 


R = RpuXp9^ X g up = 2 (K 12 + K 13 + K 23 ) = -2 In 2 A, = -2 h(T)\ 


(8.60) 


where h(T ) is the entropy of the automorphism T (8.53) 
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